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Abstract 

In this work we address the question of how oscillations are arrested in the 
mouse somitogenesis clock when the determination front reaches presomitic 
cells. Based upon available experimental evidence we hypothesize that the 
mechanism underlying such a phenomenon involves the interaction between 
a limit cycle (originated by a delayed negative feedback loop) and a bistable 
switch (originated by a positive feedback loop). With this hypothesis in mind 
we construct the simplest possible model comprising both negative and pos- 
itive feedback loops and show that (with a suitable choice of paremeters): 
1) it can show an oscillatory behavior, 2) oscillations are arrested via an 
infinite-period bifurcation whenever the different gene-expression regulator- 
inputs act together in an additive rather than in a multiplicative fashion, 
and 3) this mechanism for oscillation arrest is compatible whit plentiful ex- 
perimental observations. 

Keywords: Mathematical model, Delay differential equation, Nonlinear 
dynamics, Infinite period bifurcation 



1. Introduction 

As elegantly reviewed by Campanelli and Gedeon pQ , somitogenesis is the 
process by which vertebrate embryos develop somites, which are transient, 
repeated blocks of cells arising from the presomitic mesoderm (PSM) that 
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differentiate further into vertebrae, ribs, musculature, and dorsal dermis. 
The tail bud is a proliferative zone at the posterior end of the embryo where 
immature cells are continually added to the posterior-most PSM. As the tail 
bud grows away posteriorly, the oldest cells in the anterior PSM segment in 
groups to form lateral pairs of somites along the midline. The process stops 
when the anterior formation of somites has progressed posteriorly across the 
entire PSM (reaching the arresting growth in the tail bud) [2HI]. 

The first model attempting to explain somitogenesis was due to Cooke 
and Zeeman [5]. They postulated that the susceptibility of cells in the PSM 
to form somites oscillates between susceptible and unsusceptible, while a 
determination wavefront sweeps posteriorly across the PSM. The passing 
wavefront triggers cells to form somites, but does so only when cells are 
susceptible, i.e., when their clocks are in the correct phase of oscillation. Since 
adjacent cells are in phase, cohorts of cells are recruited in succession to form 
somites. Initially, the clock was thought to be closely linked to the cell cycle 
[B] . However, Palmerim et al. [7] discovered a gene with oscillatory expression 
in the PSM of the chick embryo, providing an alternative candidate for the 
clock. Experimental work has since identified multiple oscillatory genes in 
each of several model organisms, including mouse [8] and zebrafish [9]. 

In all of these organisms, the oscillatory gene expression in individual 
cells is coordinated throughout the PSM in order to produce spatiotemporal 
waves of mRNA and protein expression [Ml [EE], which we call the clock- wave 
[7J. Synchronized, periodic expression is observed in the tail bud with a 
frequency that matches the anterior formation of somites [HE]. Broad waves 
of expression repeatedly initiate in the posterior-most PSM and narrow while 
traveling anteriorly [31 SI [7] . The waves slow considerably as they reach the 
region of forming somites. Successive waves arriving at the anterior-most 
PSM help sequentially establish stable bands of high- low gene expression in 
several additional genes, indicating nascent somite boundaries and polarity 
[» 

Separate experiments have identified biochemical candidates for the wave- 
front [T2TfTTj . These bio- molecules exhibit graded concentration profiles across 
the PSM that shift posteriorly in synchrony with tail bud growth. A changing 
gradient level triggers mesodermal cell differentiation and somite formation 
[U S [El HE]. We call this the gradient-wavefront. 

The precise mechanism in which the clock-wave interacts with the gradient- 
wavefront, as well as their possible interactions with intercellular signaling 
mechanisms, remains unknown [H |9l [181 HH] • Many mathematical models of 



2 



the dynamics of somitogenesis have been proposed, with reviews and com- 
parisons of several prominent models available in the literature [TBI [18H2T] . 

Previous works have suggested the presence of bistable switches arising 
from positive feedback loops (formed by mutual interactions between genes 
involved in the somitogenesis process) [T^l H7J- In a recent work [21] we 
suggested that the interaction between an oscillatory regime (originated by 
a time-delayed negative-feedback loop) and a bistable domain (originated by 
a positive feedback loop involving genes under different pathways) may be 
responsible for the arrest of oscillations in presomitic cells in such a way that 
cohorts of presomitic cells stop oscillating simultaneously. Nevertheless, such 
a model falls short to explain the experimentally-observed oscillation-period 
increase as PSM cells approach the somite- formation region [71 1281 129] . 

In this paper we develop a biologically informed, yet minimally con- 
structed, mathematical model which nonetheless is capable of explaining 
most of the observed dynamic features of oscillation arrest, and provides 
a plausible mechanism for this phenomenon; namely, an infinite period bi- 
furcation. 

2. Methods 

2.1. Model Development 

Several studies on mice have demonstrated that the expression of various 
genes under the Notch regulatory pathway oscillates in presomitic cells [TT1 
1301132] . Other modeling studies have suggest that the underlying mechanism 
is a simple negative feedback loop with relatively long time delays due to 
transcription and translation of the corresponding genes [TBI [EH1 HS1 I2S]- On 
the other hand, the expression of genes under the Wnt and FGF signaling 
pathways has also been experimentally demonstrated to oscillate [HI EI21 133^ 
135] , with the same period as that of genes under the Notch pathway. There is 
evidence that such oscillation synchronization is due to the mutual regulation 
of genes under the different regulatory pathways [3"oT43"5] . The most complete 
regulatory regulatory network involving genes known to participate in the 
somitogenesis process in the mouse is reviewed in [1]. 

Based on the above discussed facts we proposed in a previous work [21] 
a genetic circuit composed of a delayed negative feedback regulatory loop, 
capable of generating an oscillatory behavior, and a positive feedback loop, 
behaving as a bistable switch. We observed that the interaction of the switch 
and the oscillator can explain the simultaneous arrest of oscillations in cohorts 
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of cells; a necessary condition for explaining somite formation. Nevertheless, 
the model failed to explain the experimentally-observed increase of the oscil- 
lation period as presomitic cells get away from the tail bud EH [29] . Under 
the consideration that such period increase is more than just an interesting 
dynamic feature, we decided to investigate its origin in the present work. 

The simplest possible gene network comprising a negative and a positive 
feedback loops is that illustrated in Figure [T]A.. Notice that we have included 
an open-loop regulatory input, k, to account for the wavefront signal. In the 
previous paper we assumed that a negative and a positive regulatory signals 
acting upon a given gene only allow its expression when the first is absent 
and the second is present. In the present work we consider that the gene 
can be expressed when either the negative signal is absent or the positive 
signal is present. Given the above considerations, the equation governing the 
dynamics of the gene circuit in Figure [T]A. is: 



In the above equation f(x) < 1 and g(x) < 1 respectively stand for the 
negative and positive feedback loops, x T denotes variable x delayed a time 
r, a is the degradation rate for proteins x, b is the strength of the negative 
loop as compared to the positive loop, and < k < 1 is the strength of the 
wavefront signal. Equation (JIJ) has been parametrized in such a way that 
all of the x input regulatory signals add up to one when they reach their 
maximal values. This parametrization guaranties that x is a dimensionless 
variable, and that x < 1 at all times. To ensure that the sum of all regulatory 
inputs is always less than one, we have further assumed that the strength 
of both the negative and the positive loops decreases proportionally to k; 
indeed, we have supposed that both of them decrease in the same amount. 
Finally, we have assumed that functions f(x) and g(x) are Hill- type functions 
of the form: 



Observe that f(x) is a monotonously decreasing function while g(x) is a 
monotonously increasing function, in accordance with the fact that they 
stand for negative and positive regulation, respectively. 




(1) 




(2) 



(3) 
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2.2. Numerical Methods 

All the solutions and numerical analyses performed on the mathematical 
model here presented were carried out using the software xppaut and the 
numerical methods available in it [39J. 

3. Results 
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Figure 1: A. Schematic representation of the gene network studied in this work. B. 
Bifurcation diagram of the mathematical model given by Eqns. (flj)-([3]), with the model 
parameters set in such a way that an infinite period bifurcation is obtained; in particular 
t p = 0. C. Plot of the oscillation frequency vs. the bifurcation parameter k, with the 
same parameter values as in B. D. Plot of % vs. t (blue line) calculated with the same 
parameter values as in C, and with k(t) as plotted (green line). E. Same as in D, but with 
T p — 20 min. 

Of all the parameters in the model we fix n = 4, which corresponds to a 
biologically plausible level of cooperativity, and analyze the model dynamics 
under variations of the rest of them. We start by studying the behavior of 
the negative feedback loop by itself, and thus we set b = 1 and k — 0. As 
previously reported, we found an oscillatory behavior for a vast set of a, k\ 
and r n values. Moreover, the smaller the value of r n , the larger the value of a 
and the smaller the value of k%, in order for the system to present sustained 



5 



oscillations. Finally, r n is the parameter that has, by far, the largest influence 
on the oscillation period. Thus, we selected r n = 0.64 hrs = 38.4 min, a = 
10 hrs -1 , and k\ = 0.24 to reproduce the 2 hrs oscillation period reported for 
the mouse somitogenesis clock. 

A B 




Figure 2: Snapshots, taken every 30 min, of the time evolution of 12 cells aligned along 
the PSM, with a separation of 30 min between consecutive cells. The time span of the 
whole simulation begins when all the cells are in the tail bud and finishes when the last 
cell is arresting its oscillations. The simulation in A was computed with r p — (the 
height of each column is proportional to the value of x in the corresponding cell), while 
the simulation in B was carried out considering r p = 20 min. 

To analyze the interaction between the oscillator and the bistable switch 
we study the model behavior under variations of parameters k and h. Param- 
eter k plays the role of a bifurcation parameter since the oscillations of x are 
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invariably arrested as the value of k is increased beyond a given threshold. 
In many cases, oscillation arrest occurs via a Hopf-like bifurcation. That is, 
the oscillation amplitude decreases continuously as k increases, until it even- 
tually becomes zero. However, for some values of parameters b and k 2 the 
oscillation arrest occurs via an infinite period bifurcation. After thoughtfully 
browsing the parameter space we found that parameter b needs to take values 
in the range [0.2, 0.4] in order to have an infinite period bifurcation. For val- 
ues larger than 0.4, the oscillations are arrested via a Hopf-like bifurcation, 
while for b < 0.4 no oscillatory behavior is observed whatsoever. Regarding 
k 2 the range of values this parameter can attain so we observe an infinite 
period bifurcation depends on the value of b. For instance, when b = 0.2 an 
infinite period bifurcation occurs whenever k 2 € [0.46,0.49], when b = 0.3 
the range of k 2 values for which an infinite period bifurcation is observed 
is [0.41,0.52], and when b = 0.4 the condition for having an infinite period 
bifurcation is that k 2 G [0.37,0.43]. In all cases, the k bifurcation value is a 
growing function of k%. In Figure [Tj3 we present the bifurcation diagram cal- 
culated with the following b and k 2 parameter values: b = 0.3 and k 2 = 0.48; 
while in Figure [lp we plot the oscillation frequency as a function of k, for 
the same parameter values. In this last plot we can see how the frequency 
goes to zero (and so the period goes to infinity) as k reaches its bifurcation 
value. 

As the embryo grows the tail bud recedes and, while doing so, leaves 
some cells behind. After a cell leaves the tail bud, the levels of the different 
morphogenic substances (Fgf8, Wnt3a, and Retinoic Acid) start changing: 
the Fgf8 and Wnt3a levels decrease, while the Retinoic Acid level increases 
|17j . In our model we account for those changes by increasing the value 
of parameter k. In Figure [Tp we illustrate the behavior of an oscillatory 
presomitic cell as the value of parameter k increases nonlinearly with time; 
starting from zero and up to a value 4/3 times its bifurcation value. Notice 
how the oscillation period increases together with k, and how, when k sur- 
passes its bifurcation value, the system completes one further oscillation and 
then reaches a steady state. 

To better understand the process of oscillation arrest along the PSM we 
simulate, following [21], the time evolution of 12 cells aligned along the PSM, 
with a separation of 30 min between consecutive cells. If the tail bud recedes 
at constant velocity and keeps leaving cell behind in a steady fashion, the 
distance between two PSM cells is proportional to the difference of elapsed 
times since their leaving the TB. Given that the oscillation period is 2 hrs, 
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the considered cell set spans a PSM region three periods long. To account for 
the cell separation, we simply assumed in our simulations that the function 
describing how parameter k decays in time is delayed in proportion to how 
much later a given cell left the TB. In other words, if k(t) describes the time 
evolution of parameter k for the first cell, kit— (i— 1)AT) is the corresponding 
function for the z-th cell, with AT = 30 min. 

The results of the simulations described in the previous paragraph are 
plotted in Figure [2]A There, we present from top to bottom snapshots taken 
every 30 min of the state of all 12 cells. In total, these snapshots account 
for five cycles. Observe that at first all the cells oscillate synchronously. 
This happens because k ~ for all of them, simulating their stay in the 
tail bud. However, as time passes and the frontmost cells leave the tail bud, 
their oscillation period decreases until they stop oscillating and the gene 
expression remains at a high level. Another feature worth noticing is the 
fact that the oscillation arrest takes place in cohorts of four cells. That is, 
the first four cells complete three cycles before stopping, while the second 
and third four-cell sets complete four and five oscillations, respectively. The 
process described above can be more clearly visualized in the movie Movl of 
the Supplementary Material. 

The origin of the time delay r n has been attributed to the sum of the times 
taken by transcription, mRNA splicing, mRNA translocation, translation, 
post-transductional modification, diffusion of transcription factors into the 
nucleus, etc. Since all of these processes also take place in the positive 
feedback regulation loop, a time delay must be also associated to it. Let us 
denote such a time delay by r p . We investigated the influence of parameter r p 
on the system behavior and found that it is qualitatively the same whenever 
t p < 20 min. In particular, the cyclic behavior disappears as k is increased 
beyond a given value, via an infinite period bifurcation. Larger t p values 
render a more complex behavior whose analysis is beyond the scope of the 
present paper. 

We calculated the behavior of the model with r p = 12 min when k in- 
creases nonlinearly with time; starting from zero and up to 0.2. This simu- 
lates a single presomitic cell from the time it is in the tail bud, until it stops 
oscillating. The results are plotted In Figure [1)5. Observe that the results 
are qualitatively equivalent to those in Figure [Tp. The most noticeable dif- 
ferences are that the oscillation amplitude is smaller when r p > 0, and that 
the gene expression level increases after the oscillations are arrested. 

To test the influence of the time delay t p on the process of oscillation 
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arrest along the PSM, we repeated the simulations in Figure [2j\, but now 
with t p = 12 min. The results are plotted in Figure |2]B. A comparison of 
Figures |2|\ and |2j3 reveals that the behavior of both models is qualitatively 
equivalent. The only noticeable differences are that, in the model accounting 
for a time delay in the positive feedback loop, the oscillation amplitude is 
smaller, the increase of the gene expression level after the oscillations are 
arrested is more notorious, and the cell cohorts are more clearly defined 
when they stop cycling. The process described above can be more clearly 
visualized in the movie Mov2 of the Supplementary Material. 

4. Conclusions 

Under the assumption that oscillation arrest in the somitigenesis clock 
of the mouse takes place due to the interaction between the oscillator and 
a bistable switch, we have constructed and analyzed the simplest possible 
model that accounts for this interaction. Although the resulting model is so 
simple that no quantitative predictions can be expected, we are convinced 
that it captures enough detail of the interactions between some of the genes 
involved in the somitogenesis processes to render general qualitative predic- 
tions. 

The results of the model agree with the following experimental observa- 
tions regarding oscillation arrest in presomitic cells: 

1. All cells oscillate in synchrony while they are in the tail bud, but their 
oscillation period increases with time after they leave the tail bud and 
enter the presomitic mesoderm. 

2. PSM cells stop oscillating in equally sized cohorts. These cohorts in 
turn give rise to somites. 

3. The expression levels of the genes under the Notch pathway increases 
after they stop oscillating, and they reach a stationary state of high 
expression [121 1521 130] . 

To our consideration, the above enumerated results validate the model to 
a level that makes feasible the derivation of biological conclusions out from it. 
In particular, our model results suggest that oscillation arrest takes place via 
an infinite period bifurcation, that a positive and a negative feedback loop 
need to act together (in an additive, rather than in a multiplicative fashion) 
upon a single gene in order for this bifurcation to occur. The model predic- 
tions regarding the existence of interacting positive and negative feedback 
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loops, and that regarding the additivity of both loops can be experimentally 
tested in principle. 
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